home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
QRZ! Ham Radio 8
/
QRZ Ham Radio Callsign Database - Volume 8.iso
/
pc
/
files
/
dsp
/
srffttar.z
/
srffttar
/
SRFFT
/
fftg.c
< prev
next >
Wrap
C/C++ Source or Header
|
1991-08-26
|
9KB
|
302 lines
/* ------------------------------ fftg.c ------------------------------------ */
/* */
/* Author: Eyal Lebedinsky */
/* Date: May 1990 */
/* Version: 9 June 1991 */
/* */
/* Program generator for integer DFFT. This program was developed for a very */
/* specific need and no attempt was made to make it general. Still, the basic */
/* algorithm is a standard one from the literature. For maximum performance */
/* of the integer calculation a special mechanism was employed to track shift */
/* operations that could be delayed. Coefficients are kept as 16 bits with 1 */
/* sign bit, 1 integer bit and 14 fractional bits. So a typicall multiply */
/* will do 16bit*16bits=32bits, then shift down one bit and use the HIGH 16 */
/* bits as the result. Before adding two numbers they must be scaled down one */
/* bit to avoid overflow. Adding more numbers calls for more scaling. The */
/* delayed scaling allows the intermediate array to hold numbers that are */
/* scaled to a varying degree until a clash occurs. Well, this space is too */
/* short to elaborate any further and if you did similar kind of work then */
/* you should know what I am talking about. */
/* */
/* usage: fftgc File Size EpName */
/* File - output file name. */
/* Size - sample order: 8 for 256 points, 10 for 1024... */
/* EpName - name of the generated function ('_' added) */
/* Now assemble the output file and then use as: */
/* */
/* #define M 8 */
/* #define N (1 << M) */
/* short x[N+1]; [input/output array] */
/* short qf[(N/2)+1]; [power spectrum] */
/* extern void far fft () */
/* */
/* [then call it as: fft();] */
/* */
/* This is a modified version of Sorensen et al, IEEE ASSP-35 no 6. */
/* */
/* - uses 3+3 complex mult. */
/* - no pre scrambling. descrambling done when power spectrum calculated. */
/* - moved to 0-origin as proper for c. */
/* - trig values in table. */
/* - the case for sin==cos treated separately. */
/* - uses delayed scaling down. */
/* - does length 256 fft with 642 mults. */
/* */
/* */
/* This program is released into the public domain. */
/* */
/*----------------------------------------------------------------------------*/
#include <stdio.h>
#include <math.h>
#define SC15 32768 /* scaling factor */
static short mm[1025]; /* fp(16, 0) */
static short ad[1025]; /* fp(16, 0) */
static int cntx[9];
static short x[1025]; /* fp(16,0) */
static short qf[513]; /* fp(16,0) */
static char *file_name, *ep_name;
extern void start_fft (), end_fft ();
extern void fft1 (), fft2 (), fft3 (), fft4 (), fft5 (), fft7 (), fft8 ();
static int mad; /* maximum adjust at this level */
static void
smad (i, j) /* set proper adjust */
int i, j;
{
while (ad[i] < (mad-j)) {
fprintf (stderr, "x[%u] >>= %d (%d-%d)\n", i, 1, mad, j);
fft1 (mm[i]);
++ad[i];
cntx[0] += 1;
cntx[1] += 1;
cntx[2] += 1;
cntx[5] += 1;
}
}
main (argc, argv)
int argc;
char *argv[];
{
int n, m, j, i, k, is, id, n1, n2, n4, n8, ind,
i0, i1, i2, i3, i4, i5, i6, i7, i8,
cc1, ss1, cc3, ss3;
float e, a, a3;
if (argc < 4) {
printf ("usage: fftg File PowerOf2 EpName [stderrFile]\n");
exit (0);
}
file_name = argv[1];
sscanf (argv[2], "%u", &m);
ep_name = argv[3];
if (argc > 4)
freopen (argv[4], "wt", stderr);
start_fft (file_name, ep_name);
/* ---------- Initialise ------------------------------------------------ */
n = 1 << m;
for (i = 0; i < 9; ++i)
cntx[i] = 0;
/* ---------- Digit reverse counter -------------------------------------- */
for (i = 1; i <= n; ++i) { /* '-1' for shifting to 0-origin */
mm[i] = i-1;
ad[i] = 0;
}
j = 1;
n1 = n - 1;
for (i = 1; i <= n1; ++i) {
if (i < j) {
mm[i] = j-1; /* '-1' for shifting to 0-origin */
mm[j] = i-1; /* '-1' for shifting to 0-origin */
}
k = n / 2;
while (k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
}
/* ---------- Length two butterflies ------------------------------------- */
mad = 0;
is = 1;
id = 4;
do {
for (i0 = is; i0 <= n; i0 += id) {
i1 = i0 + 1;
smad (i0, 0);
smad (i1, 0);
fft2 (mm[i0], mm[i1]);
++ad[i0];
++ad[i1];
cntx[0] += 2;
cntx[1] += 2;
cntx[2] += 2;
cntx[3] += 2;
}
is = 2 * id - 1;
id = 4 * id;
} while (is < n);
/* ---------- L shaped butterflies --------------------------------------- */
n2 = 2;
for (k = 2; k <= m; ++k) {
++mad;
n2 = n2 * 2;
n4 = n2 / 4;
n8 = n2 / 8;
e = 6.2831853071719586 / n2;
is = 0;
id = n2 * 2;
do {
cc1 = (SC15/2) * sqrt (0.5);
for (i = is; i <= n - 1; i += id) {
i1 = i + 1;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
smad (i1, 0);
smad (i3, 1);
smad (i4, 1);
fft3 (mm[i1], mm[i3], mm[i4]);
ad[i1] += 1;
ad[i3] += 2;
ad[i4] += 2;
cntx[0] += 3;
cntx[1] += 3;
cntx[2] += 5;
cntx[3] += 4;
if (n4 != 1) {
i1 = i1 + n8;
i2 = i2 + n8;
i3 = i3 + n8;
i4 = i4 + n8;
smad (i1, 1);
smad (i2, 0);
smad (i3, 1);
smad (i4, 1);
fft4 (mm[i1], mm[i2], mm[i3], mm[i4], cc1);
ad[i1] += 2;
ad[i2] += 1;
ad[i3] += 2;
ad[i4] += 2;
cntx[0] += 4;
cntx[1] += 4;
cntx[2] += 2;
cntx[3] += 6;
cntx[4] += 2;
}
}
is = 2 * id - n2;
id = 4 * id;
} while (is < n);
a = e;
for (j = 2; j <= n8; ++j) {
a3 = 3 * a;
cc1 = (SC15/2) * cos (a);
ss1 = (SC15/2) * sin (a);
cc3 = (SC15/2) * cos (a3);
ss3 = (SC15/2) * sin (a3);
/* a = j * e;*/
is = 0;
id = 2 * n2;
do {
for (i = is; i <= n - 1; i += id) {
i1 = i + j;
i2 = i1 + n4;
i3 = i2 + n4;
i4 = i3 + n4;
i5 = i + n4 - j + 2;
i6 = i5 + n4;
i7 = i6 + n4;
i8 = i7 + n4;
smad (i1, 0);
smad (i2, 0);
smad (i3, 2);
smad (i4, 2);
smad (i5, 0);
smad (i6, 0);
smad (i7, 1);
smad (i8, 1);
if (ind = (ad[i3] < (mad-1))) {
++ad[i3];
++ad[i4];
cntx[2] += 1;
}
fft5 (mm[i1], mm[i2], mm[i3], mm[i4],
mm[i5], mm[i6], mm[i7], mm[i8],
ss1-cc1, -ss1-cc1, cc1,
ss3-cc3, -ss3-cc3, cc3, ind);
ad[i1] += 1;
ad[i2] += 1;
ad[i3] += 2;
ad[i4] += 2;
ad[i5] += 1;
ad[i6] += 1;
ad[i7] += 2;
ad[i8] += 2;
cntx[0] += 8;
cntx[1] += 8;
cntx[2] += 4;
cntx[3] += 18;
cntx[4] += 6;
}
is = 2 * id - n2;
id = 4 * id;
} while (is < n);
a += e;
}
}
++mad;
for (i = 1; i <= n; ++i)
smad (i, 0);
fft7 (0, mm[1]);
for (i = 1; i < n/2; ++i)
fft8 (i, mm[i+1], mm[n-i+1]);
fft7 ((n/2), mm[n/2+1]);
end_fft (file_name, ep_name);
for (i = 0; i < 9; ++i)
fprintf (stderr, " %5u", i);
fprintf (stderr, "\n Load Store Shift Add Mult Adj\n");
for (i = 0; i < 9; ++i)
fprintf (stderr, " %5u", cntx[i]);
fprintf (stderr, "\n");
exit (0);
}